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A manifestation of a novel type of hysteresis related to the parametric resonance in the system 
of coupled Josephson junctions is demonstrated. Opposite to McCumber and Steward hysteresis, 
we find that the width of this hysteresis is inversely proportional to the McCumber parameter and 
depends also on coupling between junctions and the boundary conditions. An investigation of time 
dependence of the electric charge in superconducting layers allow us to explain the origin of this 
hysteresis by different charge dynamics for increasing and decreasing bias current processes. The 
effect of wavelength of the longitudinal plasma waves created at the resonance on the charging of 
superconducting layers is demonstrated. We found a strong effect of the dissipation in the system 
on the amplitude of the charge oscillations at the resonance. 



The hysteresis features of single Josephson junction 
(JJ) were studied by McCumber and Steward long time 
ago. m Particularly, it was shown that the width of 
the hysteresis is directly proportional to the McCumber 
parameter Pc related to the dissipation parameter {3 by 
/3c = In compare with the single Josephson junc- 

tion, the system of the coupled Josephson junctions has 
a multiple branch structure and the definition of the re- 
turn current is more general now: the system can re- 
turn to the zeroth voltage state from any branch. The 
branches have a breakpoint (BP) and a breakpoint region 
(BPR) before transition to another branch. [3[ The BP 
current characterizes the resonance point at which the 
longitudinal plasma wave (LPW) with a definite wave 
number k is created by Josephson oscillations in stacks. 
Hysteresis features of intrinsic JJ in high temperature 
superconductors have a wide interest also due to the ob- 
served powerful coherent radiation from such system. [4] 
In Ref. [H^ the authors summarized the experimental re- 
sults and stressed that the strong emission was observed 
near the unstable point of the return current in the uni- 
form voltage branch. The radiation is related to the same 
region in the current voltage characteristics (CVC) where 
the BP and the BPR were observed and it made the phase 
dynamics investigation of intrinsic JJ corresponding to 
these parts of CVC an actual problem today. 

Since thickness of superconducting layer (S-layer) in 
the intrinsic JJ is comparable to the Debye screening 
length, the S-layers are in the nonstationary nonequi- 
librium state due to the injection of quasiparticle and 
Cooper pairs. The charge neutrality in S-layers 

is locally broken and this charging effect modifies the 
Josephson relation between voltage and the phase differ- 
ence. The question concerning the value of the electric 
charge in S-layer and its maximum realized at the para- 
metric resonance is not investigated yet. 

In this paper we study the phase dynamics in coupled 
Josephson junctions. A novel type of hysteresis related to 
the parametric resonance in this system is demonstrated. 
We show that the width of this hysteresis is inversely pro- 



portional to the McCumber parameter and depends on 
coupling parameter of the system and the boundary con- 
ditions. The origin of this hysteresis is related to the 
different charge dynamics for increasing and decreasing 
bias current processes. We discuss the question concern- 
ing the maximal electric charge in S-layers realized at the 
resonance and show that it depends on the relation be- 
tween the wavelength of LPW and the period of lattice. 
We demonstrate a strong effect of the dissipation in the 
system on the coefficient of the exponential growth of the 
maximal electric charge in S-layers. 

First we discuss the phase dynamics in the system 
of coupled J J. The CVC of JJs are numerically calcu- 
lated in the framework of capacitively coupled Josephson 
junction model with diffusion current (CCJJ+DC) 0. 
The system of equations for the gauge- invariant phase 
differences ipi{T) = 6>^+i(r) - Oi{t) - ^ J^^^ dzAz{z,r) 

,2 

between S'-layers in this model has the form ^(fi = 
{l-a\/'^^^){I-sm{(pi)-/3^) where Oi is the phase of the 
order parameter in the S-layer /, is the vector poten- 
tial in the barrier, a and /3 are coupHng and dissipation 
parameters, respectively, and V^^V^ = fi+i + fi-i - "^fi- 
In our simulations we use both periodic and nonperiodic 
boundary conditions (BC). At nonperiodic BC it is sug- 
gested that the first and the last S-layers are in contact 
with normal metals and their effective width sq and sn 
may be extended due to the proximity effect into attached 
metals. Nonperiodic BC are characterized by parameter 
7 = s/sq = s/sn and the equations for the first and 
last layers in the system for phase differences are differ- 
ent from the equations for the middle S-layers 0, HI- We 
solve the system of dynamical equations for phase differ- 
ences using the fourth order Rungge-Kutta method. We 
use a dimensionless time r = tujp^ where ujp is the plasma 
frequency ujp = ^j2elcJhC ^ Ic is the critical current and 
C is the capacitance. In our simulations we measure the 
voltage in units of Vq = hup/ (2e) and the current in units 
of the /c. The details concerning numerical procedure are 
given in Ref. 9. 
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FIG. 1: Demonstration of the hysteresis behavior in the para- 
metric resonance region in the outermost branch of the CVC 
for the stack with 9 JJ at a = 1, /3 = 0.2 and 7 = 0.5 . The 
arrows show the direction of the bias current sweeping. The 
inset (a) shows the total branch structure in the CVC of this 
stack and BP location. The inset (b) shows a change of the 
hysteresis width /s* — Ib with dissipation parameter (3 and 
coupling parameter a. The inset (c) shows the same with 
variation of the nonperiodic parameter 7. 

The results of simulation of CVC and its features near 
parametric resonance region are presented in Fig. [TJ The 
inset (a) shows the total branch structure of the CVC for 
the stack with 9 JJ at a = 1, /3 = 0.2 and 7 = 0.5. The 
CVC shows the following features: (i) a jump at I/Ic = 
1.0 from the zero voltage branch to the outermost branch 
with all junctions in the rotating state; (ii) practically 
linear dependence of the voltage on bias current at / > 
Ic] (iii) multiple branching in the hysteresis region. The 
circle and arrow with letter B show the BP location on 
the outermost branch. 

Fig. [1] demonstrates a hysteresis in the outermost 
branch of the CVC. It is obtained by decreasing the bias 
current till some point in the BPR (curve 1), then we 
increase the current to pass the resonance region (curve 
2). The arrows show direction of the bias current chang- 
ing. The hysteresis is characterized by its width Ib^—Ib^ 
where Ib is the value of breakpoint current in decreasing 
current process and J^* is a characteristic current value 
in the increasing current process. 

The dependence of the hysteresis width on the dis- 
sipation parameter (3 is shown in the inset (b). The 
width is increasing with parameter P (i.e., it is decreas- 
ing with the McCumber parameter). As we mentioned 
above, the McCumber and Steward hysteresis demon- 
strates an opposite behavior for single J J. They obtained 
that the return current Ir (which characterizes the value 
of hysteresis) decreases with increasing of the McCum- 
ber parameter. [iq| So, we have observed a novel type of 
hysteresis related to the parametric resonance in coupled 
J J. The insets (b) and (s) show as well an increase of the 



hysteresis width Ib^ — Ib with the coupling parameter a 
and parameter of nonperiodicity 7, respectively. 

To make clear the origin of this hysteresis, we study 
time dependence of the charge oscillations in the S-layers. 
Using Maxwell equation div{eeoE) =Q, where e and eo 
are relative dielectric and electric constants, we express 
the charge density Qi = Qoa{Vi-^i — Vi) in the S-layer 
/ by the voltages Vi and Vi-^i in the neighbor insulating 
layers, where Qo = esoVo/r'j)^ and ro '^^ Debay screening 
length. The charge dynamics in the S-layers determines 
the features of current voltage characteristics of the cou- 
pled Josephson junctions. 
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FIG. 2: Difference in the charge-time dependence and CVC 
in (a) decreasing current process and (b) increasing current 
process. The thick curves show the CVC. 

Solution of the system of dynamical equations for 
phase differences gives us the voltages as a functions 
of time Vi{t) in all junctions in the stack, and it al- 
lows us to investigate the time dependence of the charge 
in each S-layer. Here we investigate the charge-time 
dependence for two processes: decreasing (Figs. 2(a), 
3, 4, 5) and increasing (Fig. EJb)) the bias current 
through the stack. The recorded time is calculated as 
tr = tUp + Tjn {Iq — foi" decreasing bias current pro- 

cess. For increasing current process (Fig.[2jb)) we record 
the time dependence at bias current value / during the 
time interval {tr^tr — Tm)- We put mostly Tm = 1000, 
Sr = 0.05 and SI = 0.0001 in our simulations. 

In Fig. [2] the time dependence of the charge in S-layer 
of the stack with 9 coupled JJata = l, /3 = 0.2 and 
periodic BC is combined with the CVC of the outermost 
branch. Here and further we present the charge oscilla- 
tions in the first S-layer of the stack: the features of the 
charge oscillations we are interested in this papers are the 
same in the all other layers. Fig. [2fa) shows the charge- 
time dependence when the current approaches the reso- 
nance point in decreasing current process, while Fig.[2fb) 
presents it when the bias current is increased. We can 
see that the charge on the S-layer in Fig. EJb) disappears 
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at different value of current in compare with the current 
value in decreasing process (Fig.[2fa)). We consider that 
it is an origin of the observed hysteresis, which is related 
to the parametric resonance in this system. 
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FIG. 3: Demonstration of the absence of fine structure in 
charge-time dependence and CVC for the stack with 10 junc- 
tions ata = l, /3 = 0.2 and periodic BC. The filled squares 
mark the bias current's steps in CVC. Inset (a) shows the 
charge distribution among the layers. Inset (b) illustrates the 
beginning of the growing region of the charge in the S-layer. 



FIG. 4: Charge-time dependence and CVC at a = 3, /3 = 0.3 
and periodic BC: (a) for a stack with 10 J J; (b) for a stack 
with 9 J J. The inset shows the charge distribution among the 
layers corresponding to the 27r/3(i-mode in the stack with 9 
JJ. 



Let us now discuss the question concerning the am- 
plitude of the electric charge oscillations in S-layer at 
parametric resonance (the parametric resonance corre- 
sponds to the BP on the outermost branch of CVC). Does 
its maximal value depend on the wavelength of created 
LPW? The Josephson oscillations excite the LPW with 
k = n/d (TT-mode, wavelength A = 2d) at the parametric 
resonance in the stack with even number of J J at a = 1 , 
/3 = 0.2 and periodic BC.[11] But in the stacks with odd 
number N of JJ the wave number depends on N and it 
is equal to k = 7t{N — l)/dN. In the stacks with even 
number of J J the resonance is " pure" , i.e., no additional 
fine structure appears in CVC.[12| So, it is interesting to 
compare the maximal value of the electric charge realized 
in S-layers in this case with the case A ^ nd^ where n is 
integer number. 

In Fig. [3] we present time dependence of the electric 
charge in the S-layer for the stack with 10 J J combined 
with the outermost branch of CVC (the corresponding 
axis are shown by arrows). We found that the bias cur- 
rent interval where the growing region (the beginning of 
that region is demonstrated in the inset (b)) of the elec- 
tric charge in S-layers is observed, is shorter now in com- 
pare with N=9 case, where the LPW with k = Sn/dd 
is created at the same values of a and p. Compare this 
figure with Fig. [21 we can see that the amplitude of the 
charge is larger for the stack with even number of J J. 
The inset (a) illustrates the charge distribution among 
the layers and confirms the 7r-mode of LPW. 

As we mentioned above, the wavelength of the LPW 
depends on the values of dissipation and coupling 



parameters. [3] So, we can compare the stacks with 10 
and 9 JJ when another LPW are created and test the 
idea concerning the maximal amplitude of charge oscilla- 
tions. In Fig. m the time dependence of the charge in the 
first S-layer ata = 3,/3 = 0.3 and periodic BC combined 
with the CVC of the outermost branch is presented. 

Fig. m^a) shows this dependence for the stack with 
10 junctions. At chosen parameters the LPW with 
k = 37r/5dis created. In case = 9 (Fig.lU^b)) the inset 
ihustrates the charge distribution among the layers corre- 
sponding to the (27r/3(i)-mode (A = 3d). Also we can see 
that the charge value on the S-layers in k = 2T:/3d case 
is larger than in case oi k = 37r/5(i, which is the same 
result as we got before. In addition to that, we tested 
the cases for \ = Ad and A = 5(i ( not presented here) 
and they supported our idea. So, we may conclude that 
at fixed a and {3 the charge value in S-layers is larger for 
the stacks with ^'pure^' parametric resonance where the 
LPW with X = nd is created. 

To demonstrate the character of the charge amplitude 
increasing in the growing region, we enlarged in Fig. [5] 
the charge-time dependence for the stack with 10 J J at 
a = 1 and /3 = 0.2. In inset (a) we present the time 
dependence of the amplitude of Q^/Qo in the logarith- 
mic scale. The values of amplitude were taken arbitrary 
at some time moments in total growing region (examples 
are shown by circles). We found two parts in this depen- 
dence: exponential part and transition part (marked by 
double arrow) before jump to another branch. In tran- 
sition part the amplitude demonstrates a sharp increase 
in short time interval in compare with exponential part. 
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FIG. 5: Demonstration of the transition part (shown by dou- 
ble arrow) in charge-time dependence for the stack with 10 
J J. Inset (a) shows the amphtude of charge oscillations in the 
logarithmic scale. Insets (b) and (c) demonstrates the results 
of FFT analysis of voltage V{t) and charge Qi{t)/Qo time 
dependence in the exponential growth part. 
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FIG. 6: Influence of dissipation magnitude on time depen- 
dence of the charge oscillation amplitude in the logarithmic 
scale for the stack with 10 junctions at a = 1. 



efficient K. The value of K is determined by the wave 
number of LPW created at resonance. For all investi- 
gated stacks with even number of JJ (in our simulations 
we checked the stacks with N in the interval (6,14)) at 
a = 1, /3 = 0.2 and periodic BC we obtained the same 
value of = 0.001. 

In the supplement to this letter we include the avi files 
which present the animation of the charge dynamics in 
the different time moments in the growing region. These 
animation files demonstrate the standing tt — mode of 
created LPW (the charge on the nearest neighbor layers 
has the same value and opposite sign) in the exponential 
part and its modification in the transition region. 

As a summary, a manifestation of a novel type of hys- 
teresis related to the parametric resonance in the system 
of coupled Josephson junctions is demonstrated. The 
width of this hysteresis is inversely proportional to the 
McCumber parameter and depends on coupling between 
junctions and the boundary conditions. The origin of 
this hysteresis is related to the different charge dynamics 
for increasing and decreasing bias current processes. We 
consider that these features are common for the systems 
demonstrating the parametric resonance. These features 
can be used to develop new methods for determination of 
coupling and dissipation parameters of the system. We 
show that the maximal value of electric charge ampli- 
tude realized in superconducting layers at the resonance 
depends on the wavelength of the created LPW. A strong 
effect of the dissipation in the system on the width of the 
parametric resonance is demonstrated. 

We thank I. Rahmonov, M. Hamdipour, H. El Sam- 
man, S. Maize, M. Elhofy and Kh.Hegab for helpful dis- 
cussions and support this work. 



Insets (b) and (c) show results of FFT analysis of voltage 
V{t) and charge Qi{t)/Qo time dependence in the expo- 
nential grows part. They prove the parametric resonance 
condition ujj = 2ujlpw- In transition part this condition 
is broken. Writing the expression for electric charge by 
Qi/Qo = exp(Kt), we find K = 0.001. 

The Fig. [6] illustrates the influence of the dissipation 
magnitude {/3 = 0.2, 0.1, and 0.05) on time dependence of 
the charge oscillation amplitude in the logarithmic scale 
for the stack with 10 junctions at a = 1. From this 
figure we note the following features which are observed 
with decrease in P (increase the McCumber parameter): 
i) the growing region is getting shorter; ii) the width of 
transition part is decreasing; iii) the coefficient of the 
exponential growth K is increased. We come to the im- 
portant conclusion that the parametric resonance features 
depend strongly on dissipation in the system. The width 
of the growing region is inversely proportional to the co- 
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